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Abstract 

We evaluate the virial coefficients Bk for fc < 10 for hard spheres in dimensions D = 2, • • • , 8. Virial 
coefficients with k even are found to be negative when 7? > 5. This provides strong evidence that the 
leading singularity for the virial series lies away from the positive real axis when D > 5. Further analysis 
provides evidence that negative virial coefficients will be seen for some A; > 10 for D = A, and there is a 
distinct possibility that negative virial coefficients will also eventually occur for D = 3. 

Keywords: hard spheres, virial expansion. 



1 Introduction 

The hard sphere gas of particles of diameter a \n D dimensions defined by the two body potential 

is one of the oldest and most studied systems in statistical mechanics. Nevertheless after more than a century 
after the initial work of van der Waals PP , Boltzmann |2] , and van Laar |3] who computed the coefficients up 
through ,84 in three dimensions for the low density virial expansion of the pressure 

p °c 

Y^^Y. ^kp" with B, = 1 (2) 
^ fc=i 

there are still basic features of this system which are unresolved and controversial. Chief among these properties 
are 1) the radius of convergence of the virial series, 2) the question of the occurrence of negative virial 
coefficients in 2 and 3 dimensions, and 3) the analytic demonstration of the freezing phase transition which 
was first seen in computer experiments ^ E] in the late 1950s. 

The analytic information available for the hard sphere gas is extremely sparse. Indeed, outside of B2 
and B3 which are elementary to compute and the original computation of -B4 in D = 3 by Boltzmann and 
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van Laar, the only other analytic computations of virial coefficients are of B4 in D = 2 done simultaneously 
by Rowlinson 6 and Hemmer (7j in 1964, and the recent computation by the present authors [5| of B4 for 
D = 4, 6, 8, 10, 12, and by Lyberg for odd dimensions £> = 5, 7, 9, 11. 

All other computations for the hard sphere gas are by means of computer. The fifth virial coefficient 
was first calculated numerically in the 1950s for hard discs by Metropolis et al. JOI and for hard spheres 
by Rosenbluth and Rosenbluth ^J. An extensive and systematic effort to numerically calculate the virial 
coefficients B5, Bq, and B^ for two and three dimensions was carried out by Ree and Hoover [T^ 1 1 31 [H] during 
the 1960s. The coefficient B^ hi D — 2 and 3 was computed by Janse van Rensburg in 1993. Vlasov, You, 
and Masters ^H] recently recalculated B^ and Bg, and Labik, Kolafa, and Malijevsky |17[ll8j have recalculated 
the virial coefficients to B^ and also calculated Bq. All these virial coefficients are positive. 

The study of virial coefficients for dimensions greater than 3 was initiated in 1964 by Ree and Hoover ^H] 
who computed B4 for Z) = 4, • • • , 11 and found for D > 8 that B^ is negative. This was the first time that 
any negative virial coefficient had been seen for the hard sphere system. The coefficients B^ and Bq for = 4 
and 5 were computed by Bishop, Masters, and Clarke ^0 in 1999, and Bishop, Masters, and Vlasov 7^ have 
recently calculated By in dimensions four and five, and Bg, in four dimensions. The present authors recently 
extended the study of B4, B^, and Bq up to dimensions D = 50. In this study it was found that B^ is 
always positive but the Bq is negative for D > 6. 

There are two very important observations to be made about these computations. 

The first observation is that the virial coefficients B2 through Bg in D — 2,3 have been analyzed by 
many authors in an attempt to find an approximate equation of state. These approximations are reviewed 
in Subsection 14.31 Without exception all of these studies have the remarkable feature that they show no 
singularity at the density at which the computer studies indicates that the system freezes. Indeed some of 
the approximate equations of state are analytic for real positive densities greater than the density of closest 
packed discs and spheres. One interpretation of this is that the first eight virial coefficients are not sufficient 
to see the true asymptotic behavior of Bk for large k. 

The second observation is that because B4 changes sign between D = 7 and D — 8, and Bq changes sign 
between D = 5 and 6 it may be for large k that Bk can become negative for dimensions smaller than 5. In 
particular if for D = 2 or D = 3 there were a value of k such that Bk changed sign then the approximate 
equations of state obtained from the first eight virial coefficients would be wholly inadequate to obtain the 
radius of convergence of the virial series. 

In this paper we address the questions of the radius of convergence and and the negativity of the virial 
coefficients by numerically computing Bj for D = 6,7,8, Bg, for D — 5,--- ,8, and Bg and Biq for D = 
2, • • • ,8. Our results are given in Tabled For the Monte-Carlo integration we use the formulation of Ree 
and Hoover I12L I13| . For hard spheres it is well known that in any dimension D there may be some Ree- 
Hoover diagrams which vanish identically for geometric reasons. The determination of the number of such 
geometrically excluded diagrams is an interesting problem in its own right and our results for this are given 
in Table El 

In Section 13 we discuss the virial expansion and the theoretical framework for the problem of calculating 
virial coefficients for hard spheres. In Section|21we explain our method of computation with particular emphasis 
on what needs to be done to handle the 4,980,756 Ree-Hoover diagrams contributing to Bio in an efficient 
fashion. In Section 0] we provide information about the hard sphere system relevant to the understanding of 
the asymptotic behavior of the virial series. In Sectional we address the question of the large k behavior of Bk- 
We propose in Subsection 15.11 two complementary criteria which Bk should satisfy in order to be considered 
asymptotic. We discuss the asymptotic behavior of the virial series in 15.21 from a ratio analysis. In Subsection 
15.31 we apply the methods of Fade approximants and differential approximants to the ten known coefficients. 
We conclude in Section with a summary of key results. 

The results of this work as well as the papers of Clisby and McCoy [HI 122] are included in the dissertation 
of Clisby [21. 
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D 




B4/B2 






B7/B« 


Bs/Bl 




B10/B2 


2 


0.782004- •• 


0.53223180 •• • 


0.33355604(1)* 


0.1988425(42) 


0.1148728(43) 


0.0649930(34) 


0.0362193(35) 


0.0199537(80) 


3 


0.625 


0.2869495 • • • 


0.110252(1)' 


0.03888198(91) 


0.01302354(91) 


0.0041832(11) 


0.0013094(13) 


0.0004035(15) 


4 


0.506340- • • 


0.15184606 •• • 


0.0357041(17) 


0.0077359(16) 


0.0014303(19) 


0.0002888(18) 


0.0000441(22) 


0.0000113(31) 


5 


0.414063 - •• 


0.0759724807- • • 


0.0129551(13) 


0.0009815(14) 


0.0004162(19) 


-0.0001120(20) 


0.0000747(26) 


-0.0000492(48) 


6 


0.340941 • • • 


0.03336314 •• • 


0.0075231(11) 


-0.0017385(13) 


0.0013066(18) 


-0.0008950(30) 


0.0006673(45) 


-0.000525(16) 


7 


0.282227- • • 


0.00986494662 • • • 


0.0070724(10) 


-0.0035121(11) 


0.0025386(16) 


-0.0019937(28) 


0.0016869(41) 


-0.001514(14) 


8 


0.234614- • • 


-0.00255768- • • 


0.00743092(93) 


-0.0045164(11) 


0.0034149(15) 


-0.0028624(26) 


0.0025969(38) 


-0.002511(13) 



Table 1: Numerical values of virial coefficients. Values for Bi D > 5, Bg, D > A, Bg, and Bio are new, and other values improve on published 
literature results for B^ and higher except for the results for B^ for £> = 2, 3 which are due to Kratky j24j . 



Table 2: Number of Mayer and Ree- Hoover integrals 



Order 





2 


3 


4 


5 


6 


7 


8 


9 


10 


Mayer 


1 


1 


3 


10 


56 


468 


7123 


194066 


9743542 


RH 


1 


1 


2 


5 


23 


171 


2606 


81564 


4980756 


RH/Mayer 


1 


1 


0.667 


0.500 


0.410 


0.365 


0.366 


0.420 


0.511 


RH, D = 1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


RH, D = 2 


1 


1 


2 


4 


15 


73 


>647 


>8417 


>110529 


RH, D = 3 


1 


1 


2 


5 


22 


161 


>2334 


>60902 




RH, D = 4 


1 


1 


2 


5 


23 


169 


>2556 


>76318 





2 Virial Expansion 

To calculate virial coefficients we use the reformulation of the Mayer series |23 carried out by Ree and 
Hoover 

In the Mayer formulation virial coefficients Bf, are given as the sum of integrals which may be represented 
by biconnected graphs of k points. We define a k point graph G as a set of vertices V = {vi, i = 1, • • • ,k} 
together with a set of edges E whose elements are pairs {u, v) where u, veV. We will only be interested in 
undirected graphs, and not in generalizations such as directed graphs or graphs with multiple edges and loops. 
A graph is biconnected if there are no vertices whose removal would result in a disconnected graph, which 
are more concisely known as articulation points. Then the virial coefficient at order k is given in the Mayer 
formalism by 

= ^ E C{G)SiG) (3) 

where S{G) is the value of the integral represented by G, is the set of all labeled, biconnected graphs, 

is the set of all unlabeled biconnected graphs. C(G) is a combinatorial factor which is defined as the total 

number of distinguishable labelings of a graph, otherwise given as: 

for a graph with k vertices, and |Aut(G)| is the cardinality of the automorphism group of the graph. 

Each graph may be identified with an integral in which the vertices represent coordinates in Z3-dimensional 
Euclidean space, and a bond between vertices i and j represents the function 

fin - rj) = exp (-[/(r; - rj)/fcBr) - 1 (5) 

For example, the integral corresponding to the graph in Fig. ^ may be written as 

5(Gex) = lim / d''rid^r2d^r3d^r4/(ri2)/(r23)/(ri4)/(r34) 

= J d^rid^r2d^r3/(ri2)/(r23)/(ri)/(r3) (6) 

where ry = ri — rj, and V is the volume to be integrated over. In the second expression for the integral, is 
defined as the origin. Naturally the value of the integral does not depend on the labeling chosen, and this is 
the reason that in Eq. Olthe sum over labeled graphs cpn be converted to a sum over unlabeled graphs. 




Figure 1: Typical labeled graph, Gcx 

A useful re-summation was performed by Ree and Hoover |12l I13| by introducing the function 

/(r) = 1 + /(r) = exp {-U{v)/kBT) (7) 

and then expanding each Mayer graph by substituting i — f — f for pairs of vertices not connected by / bonds. 
This method was previously used by Percus and Yevick [201 in comparing the exact values of the fourth and fifth 
virial coefficients with coefficients obtained from the Percus- Yevick equation, and by Percus [27) in discussing 
the derivation of the Percus- Yevick equation. One can see that upon performing this re-summation on the 
Mayer series that vertices in the resultant graphs will all be mutually connected by either / bonds or / bonds, 
and that graphs must be biconnected with respect to / bonds. In this paper, we adopt the convention that 
graph edges only correspond to / bonds; in other contexts it can be very useful to refer to the complement 
graphs in which the / bonds are edges. Now there is an additional factor for each labeled graph, called the 
star content by Ree and Hoover, which may be either positive or negative, and is equal to the total number 
of biconnected subgraphs with an even number of edges removed (including zero), subtract the number of 
biconnected subgraphs with an odd number of edges removed. Thus the Ree-Hoover expression for the virial 
coefficient Bk is 

= E SC(G)5(G) 



fc! 



Ge.i 



1 - k 
fc! 



J2 G(G)SC(G)S'(G) 



(8) 



where SC(G) is the star content and S{G) is the Ree-Hoover integral value of the graph G. As a consequence of 
the fact that for every 1 expanded there is a positive and negative coefficient, in the Ree-Hoover re-summation 
only the complete graph which has no I's contributes a net non-zero amount. Thus 



E SC(G) - 1 



(9) 



Now, to compare with the Mayer integral formalism of Eq. |3 we give the expression for the Ree-Hoover 
integral represented by Fig. 



5(Ge 



lim - / d^rirf^r2d^r3d^r4/(ri2)/(ri3)/(r23)/(ri4)/(r24)/(r34) 

V— ►oo V J 

d^rid^r2d^r3/(ri2)/(ri3)/(r23)/(ri)/(r2)/(r3) 



(10) 



In the case of hard spheres, the potential is given by Eq. ^ and so /(r) and /(r) are particularly simple: 



/(r) 



-1 


|r| 


< a 





|r| 


> cr 





|r| 


< a 


+1 


Irl 


> a 



(11) 



(12) 



The / bonds force vertices to be close together, while / bonds force vertices apart. These competing conditions 
may mean that for some graphs contributing to Bk in dimension D < k — 1 there are no point configurations 
that satisfy the conditions, and hence the corresponding integral would be zero for geometric reasons. Any 
configuration of points in _D-dimensional Euclidean space can contribute to at most one Ree-Hoover diagram, 
whereas it may contribute to many Mayer diagrams, and in particular for all k point graphs G 

\SiG)\ > \SiG)\ (13) 

The inequality is strict when one excludes the complete star graph. This leads to the primary advantage of 
the Ree-Hoover formulation over the original Mayer formulation: it greatly reduces the degree of cancellation 
between positive and negative diagrams. 



3 Method 

In Section 13.11 we give a broad description for the Monte-Carlo algorithm used to perform the numerical 
integration for our results of Table ^ In Section 13.21 we define some graph theory notation that is necessary 
for our detailed discussion of the calculation of the star content of all unlabeled graphs of up to ten points in 
Section 13.31 the determination of a minimal set of spanning trees to generate configurations for the Monte- 
Carlo algorithm in Section and the unlabeling factor for graphs of up to ten points in Section 1^51 Lastly, 
in Section ITHl we give specific details of the implementation of the algorithm. 

3.1 Monte-Carlo Algorithm 

We calculate virial coefficients via the hit or miss Montc-Carlo integration algorithm described by Ree and 
Hoover ^2). Monte-Carlo integration is appropriate for high dimensional integrals, which is certainly the case 
for the integrals we calculate as in general, the calculation of Bk for D dimensions is a (fc — 1)D dimensional 
integral. 

If we need to integrate a function F{x) over a complicated, high dimensional domain Di, the hit or miss 
Monte-Carlo method works by enlarging the domain of integration to a simple, easily characterized domain 
D2 and setting the function to zero in 1)2 \-Di. We then select points uniformly from within D2 using a pseudo 
random number generator, and the integral we are interested in is given as the volume of the enlarged domain 
1/(1)2) multiplied by the mean value of F{x) over Z?2- 

/ dx F{x) = dx F{x) Di C D2, F{x) =0 Va; e D2 \ Di 

Jdi J D2 

1 ^ 

= V{D2) \i^^-Y.F{x,) x^eD2 (14) 

The Ree-Hoover integral of a given k point biconnected graph G has a trivial integrand which is ±1, but 
the domain is all configurations (ri, r2, • ■ ■ , r^) where G and the resulting point set may be identified 
with G. To calculate the integral, we first identify a spanning tree of the graph, and use all Ree-Hoover 
diagrams which are spanned by this tree as the enlarged domain. To generate random configurations in the 
Monte-Carlo algorithm we place the first particle at the origin, and then successively place the succeeding 
particles randomly within unit balls according to the generating spanning tree. Once we have the position of 
all particles we check the inter-particle distances to determine if the distances are greater than or less than 
one, and thus map the configuration to a graph. The volume of the enlarged domain is straightforward to 
calculate, as each edge of the spanning tree contributes a factor of the volume of the unit ball in D dimensions, 
which may also be given in terms of B2 as (2^2)*^^^. 

The efficiency of the algorithm may be improved by taking in to account all labeled isomorphs of G which 
may be generated by the spanning tree, but we then must divide by this unlabeling factor to get the correct 
result. 

This method of calculating Ree-Hoover integrals may be extended to calculating the virial coefficient itself, 
which up to a factor is given by the sum of the Ree-Hoover integrals of all labeled biconnected graphs multiplied 



by their star content. Thus we require a set of spanning trees which are able to generate all biconnected graphs, 
and we also need to calculate the unlabeling factor for each of the biconnected graphs. Ree and Hoover [T^IT^ . 
Janse van Rensburg and Vlasov et al. ^| proceed by finding a set of spanning trees, and partitioning 
the biconnected graphs according to which spanning tree will be used to calculate them. This means that 
if a given spanning tree produces graphs which are not in its set then these configurations are thrown away. 
We improve on this by combining the unlabeling factors from each of the spanning trees as the sum of the 
individual unlabeling factors, and never discard a geometric graph that contributes to Bk- 

For each iteration of the Monte-Carlo calculation we use each of the set of spanning trees to generate a 
configuration, which is mapped to a graph. If the graph is biconnected then the contribution of this graph, up 
to an overall factor, is given by the star content multiplied by the number of possible labelings of the graph, 
divided by the unlabeling factor. 

In practice, for each batch we keep a tally of all of the biconnected graphs generated, and only condense 
this information using the star content and unlabeling factor values to a value for Bk at the end of the loop. 
For this work we chose a batch size of 10^ configurations, and uncertainties were calculated using the standard 
formula 



Error = 



^ (</,>-< / >)^ 
^1 - 1) 



(15) 



where there are q independent batches with value Ij. All uncertainties given in this paper are equal to one 
standard deviation, and the number of batches used in the calculation of the virial coefficients in Tableware 
given in Tabled of Appendix IXI 



3.2 Graph Theory Notation 

Here we define graph theory notation necessary for the description of the calculation of the star content and 
unlabeling factor of biconnected graphs, as well as the determination of a set of spanning trees sufficient to 
generate all biconnected graphs of a particular order. 

We will use the term subgraph of a graph G to describe graphs with the same vertex set, but with an edge 
set that is a subset of the edge set of G. We will allow G to be its own subgraph; when excluding this case we 
will refer to proper subgraphs. 

We define ^j(G) and !3§j{G) to be the set of subgraphs of G with j edges removed which are connected 
and biconnected respectively. We define the composition of S§i and SSj by applying to each element of the 
set SSjiG) and taking the union of the resulting sets. The end result is exactly ^i+j{G). Thus 

^^oSSj^SS^+^ (16) 

However, if we apply to each of the G" £ SSj (G) then we can see that each subgraph in (G) is produced 
exactly (*^"') times. Thus 

|^,(G')|- (' + ^)|^.+,(G)| (17) 

In particular, for j = 1, 

|^,+i(G)| = -^ ^ |^,(G')| (18) 

G'e3gi(G) 

We will define a minimally biconnected graph in analogy to the definition of a tree: a minimally biconnected 
graph is a graph for which no proper subgraph is biconnected. If we remove any edge from a minimally 
biconnected graph then the resulting graph will not be biconnected. 

An articulation star or clique cut-set of a graph is a proper subset of vertices that is maximally connected, 
i.e. edges connect each pair of vertices, and whose removal will result in a disconnected graph. 

Let CL be an algorithm that takes as input a graph G and defines a relabeling of the vertices of an input 
graph. Then we will call CL a canonical labeling algorithm if it can be shown that 

G'l = CL(Gi) (19) 
G^ = ^CL(G2) (20) 



G'l = G'2 ^ Gi and G2 are isomorphic 



(21) 



Thus, CL converts labeled graphs to unlabeled graphs, and this can, in the appropriate context, greatly 
reduce the scale of a problem as there are many fewer unlabeled graphs by a factor of approximately fc! at 
order k. 

This is an NP-complete problem, and as such it is believed that it is impossible to obtain an algorithm 
that can be guaranteed to determine the canonical labeling of a graph in polynomial time in the number of 
vertices k. 

An efhcient algorithm has been developed by McKay[2Hj for the purpose of calculating the automorphism 
group and also the canonical labeling of a graph. Through his own testing on random graphs he finds 
performance for the canonical labeling procedure is approximately 0{k^) with 2 < /? < 3 for different classes 
of graphs despite the fact that this is an NP-complete problem. McKay's implementation of this algorithm in 
the C programming language, nauty (no automorphisms, yes?), has been extensively used here. 



To calculate the star content of all biconnccted nine and ten point graphs, we use the formula of Ree and 
Hoover ^21 outlined in Section [2 here given in the notation of Section 



We note that \^o{G)\ = 1 if G is biconnccted, and that \^i^i{G)\ — for z > if and only if G is minimally 
biconnected. Hence we can enumerate all subgraphs with a certain number of edges of a given biconnccted 
graph, if we have this information for the subgraphs. This naturally defines a recursive procedure, where if we 
apply this enumeration algorithm to the complete star, and recursively apply it to all unlabeled subgraphs, 
we will generate this list and hence can calculate the star content. 

It is more memory efficient, however, to enumerate all subgraphs for graphs with k edges (the Ree-Hoover 
ring), then to enumerate subgraphs of graphs with fc + 1 edges, until we finally reach the complete star with 
k{k — l)/2 edges. As a by-product, all minimally biconnected graphs up to tenth order were generated and 
recorded for use in the spanning tree algorithm of Section Using this procedure the star content for all 
ten point Ree-Hoover graphs was calculated in a matter of 30 minutes. We would therefore expect that the 
star contents for Bn may be calculated in a matter of days, but the problem is significantly more demanding 
on computer hardware as a large amount of memory is needed to perform the calculation efficiently. 

We plot histograms of the number of labeled graphs versus star content for Bq, B^, Bs, Bg, and Biq, 
in Figures I^HHl of Appendix For these figures we exclude graphs with zero star content which have an 
articulation star, and it is for this reason that the value appears to be low for k — 6,7. For fc > 8 the zero star 
content value is far greater than the next largest value, and for this reason is not shown in the plots. Note 
that the histograms become sharply peaked around zero star content, and that these plots become increasingly 
symmetric as one increases the order. 

It is natural that the star contents should be peaked around zero given the constraint imposed by Eq. O 
However, it is not obvious that there should be such a strong symmetry between graphs with star contents of 
opposite sign, and it would be desirable to find a reason for this behavior. 

3.4 Spanning Trees 

To find a minimum or close to minimum set of spanning trees we proceed as follows. We take the set of 
minimally biconnected graphs obtained during the star content calculation procedure. 

We then need to enumerate all labeled spanning trees of each of the minimally biconnected graphs, and 
use nauty to help keep a list of the unlabeled spanning trees of each graph. 

The problem of enumerating all spanning trees of a graph is considerably more difficult. If a graph 
has a large number of spanning trees, then the best available algorithms take time 0{N) with sub-leading 
additive terms, where N is the number of spanning trees. There are many algorithms available which may be 



3.3 Star Content 




(22) 



8 



distinguished by the sub leading terms, or in the amount of memory required, and by whether they output 
each of the spanning trees, or output the difference between consecutive trees. 

We have implemented an algorithm [21] which uses 0{V + E + VN) time and 0{V + E) space. As we are 
only concerned with the case where V is small, the resulting increase in complexity does not matter much. 

Naturally, each biconnected graph must have at least one minimally biconnected graph as a subgraph, and 
thus if we find a set of spanning trees that can generate all minimally biconnected graphs then necessarily we 
will be able to generate all biconnected graphs. 

We could use all spanning trees as our set without any problem. The reason we do not wish to do this in 
practice is that some graphs such as the ring graph contribute significantly to the final result, but only have 
perhaps one or two of a large number of spanning trees that will be able to generate it. 

We store this information in matrix form as follows: Aij = \ then spanning tree j is a subgraph of minimally 
biconnected graph i, otherwise = 0. As posed, this is the minimum set cover problem, for which no known 
polynomial time algorithm exists. 

The minimum set cover problem may be defined as follows: given a matrix in which all entries have the 
value of either or 1, and there is at least one 1 in each row, then what are the minimum number of columns 
that can be chosen such that there is one or more I's in each row? This problem has been shown to be 
NP-complete, and the authors could not find implementations of algorithms that are efficient on average for 
random matrices in the vein of the canonical labeling algorithm of McKay. 

However, for our purposes we did not need the exact minimum, and we can apply the greedy algorithm 
which will produce a minimal set cover (in the sense that we cannot remove any spanning tree from the final 
result without leaving a graph uncovered), that will be sufficiently close to the minimum for our purposes. The 
greedy algorithm is defined in the following way: choose the column with the greatest number of f s, where if 
there is a tie we break the tie arbitrarily, then remove that column and the rows which have now been covered. 
Repeat until there are no rows remaining, at which point the algorithm terminates and returns the numbers 
of the columns which have been removed. We implemented a slight variation of this: at the beginning of the 
main loop of the algorithm, if any row has a single f then the corresponding column is chosen initially, as it 
must be included in any set cover solution. 

We apply this algorithm to find a near minimum set of spanning trees which are able to generate all 
biconnected graphs for orders fc = 5, • • • , f 0, and these are shown in Fig.[7|of Appendix lUl 

We may also optimize the procedure by choosing trees that will minimize the error in the Monte-Carlo 
procedure by efficiently calculating large diagrams. We could do this rigorously, by performing the integration 
over a range of dimensions and orders, and obtaining a different set of trees for each case. We note here that 
the Ree-Hoover ring is one of the largest diagrams to the order we calculate, and that the ring diagram and 
most other loosely packed diagrams can be efficiently generated by the tree without any leaves (it is the only 
tree required to generate as can be seen in Fig. |7|l . For this reason we choose to have multiple copies of the 
tree without any leaves in our set of spanning trees. We have compared this set with the set of all spanning 
trees and it is more efficient, but we have not performed extensive testing to justify this choice. 

We make the comment that there is no difficulty in carrying this procedure out to much higher order, as 
the number of trees and minimally biconnected graphs grow relatively slowly with order compared with the 
rate of growth of the total number of graphs. 

3.5 Unlabeling Factor 

Now that we have our minimal set of spanning trees, we need to calculate the unlabeling factor as defined by 
Ree and Hoover. 

The unlabeling factor is defined as the number of ways in which a given spanning tree may generate a 
particular labeled graph. Alternatively, we may define it as the number of labeled isomorphs of a particular 
graph which have the (labeled) spanning tree as a subgraph. 

Ree and Hoover j^^J^^, Janse van Rensburg 2SI> and Vlasov et al jE] proceed by partitioning the 
set of biconnected graphs in to those which will be calculated by spanning tree 1,2,3, etc., and determining 
the unlabeling factor for each unlabeled graph. This is a computationally intensive procedure in the same 
manner as for the star content calculation, although not quite as intractable. This is because the number of 



9 



labeled trees that are subgraphs of the complete star is fc'^~^. Thus we will need to enumerate approximately 
^fc-22fc(fc-i)/2^^l spanning trees, which grows very fast with order. 

We improve on this method in two respects; firstly, we do not partition the graphs, and allow the spanning 
trees to generate any of the biconnected graphs. In this manner we do not throw out any biconnected 
configurations that we generate. The combined effect of unlabeling factors from all of the spanning trees is 
merely the sum of the individual unlabeling factors. Secondly, in this formula we may use any unlabeled 
spanning tree more than once if we choose. As noted in Section 13.41 we use a minimal set of spanning trees 
but with repetitions of the tree with no leaves. 

We then proceed as for the star content, as Eg. 1221 applies to the operator for connected graphs "^j, with the 
difference that we need to keep track of all connected graphs (not just biconnected graphs), and we only need 
to know the unlabeling factor of each graph rather than the count of subgraphs partitioned by the number of 
edges they have. 

3.6 Details of the Monte-Carlo Procedure 

All code is written in the C programming language, and the calculation has been compiled and run on Linux 
machines with Pentium 4 processors. 

We use the random number generator of Ziff |3(J| . and in some cases we have confirmed results with the 
random number generator of Knuth |31| . 

We generate random points using the algorithm of Banerjia and Dwyer |32| for _D = 3, 4, 5, and use similar 
accept or reject algorithms for higher dimension. Alternative methods which do not require the rejection of 
any generated points, such as using Gaussian variables, may be faster for higher dimensions, but as this is not 
the bottleneck step in the algorithm we did not implement or test such alternatives. 

Given a set of A'span spanning trees which will be used to generate configurations, we wish to be able 
to generate configurations as fast as possible by minimizing the number of random points to generate. To 
generate each tree we will need to have fc — 1 random points and hence naively we need a total of (fc — l)A'span 
random points for each iteration. By using the same set of random points for each spanning tree we only need 
to generate k — 1 random points per iteration. 

The remaining detail of the Monte-Carlo calculation that needs to be discussed is how to identify graphs as 
they are being generated. The method used was to create a hash table of all canonically labeled biconnected 
graphs with non-zero star content, so that when a graph is generated one can calculate the canonical label 
and check if it is in the hash table using linear probing )31j. In practice, the canonical labeling algorithm of 
McKay [2H| is quite slow compared to the speed at which configurations are generated, and so if we possible 
we wish to avoid this step. 

For small graphs, which for the computers on which this program was run means for eight or less vertices, 
we can create an array with 2*'^'^"^)/^ elements for which the address is identified with all labeled graphs of 
k points. One can map any graph to a sequence of k(k — l)/2 bits by imposing an order on all vertex pairs, 
and then giving the bit corresponding to vertices i and j a value of 1 if (i,j) G -E, and otherwise. Thus once 
a graph is generated the canonical label can be quickly looked up in this array provided some pre-processing 
is done before the main loop of the Monte-Carlo integration. Unfortunately, as the total number of labeled 
graphs grows extremely rapidly with order, memory limitations mean that even on a supercomputer it would 
not be possible to do this for graphs with more than 9 vertices. 

For large graphs it is therefore much more difficult to perform the identification quickly, and so we attempt 
to identify if the graph has non-zero star content prior to calculating the canonical label. In particular 
contributing graphs must be biconnected, and cannot contain an articulation star, as discussed by Ree and 
Hoover To this end, a graph that has been generated is first checked to see if there are any vertices 

of degree 1, and then tested for biconnectivity, and then checked to see if it has an articulation star using 
an algorithm described below. If any of these conditions are met, then the star content is necessarily zero 
(although many graphs with zero star content do not have an articulation star); otherwise the canonical label 
is calculated using the canonical labeling algorithm of McKay ,28, . 

Ree and Hoover jJJl showed that graphs with an articulation star will have zero star content, although 
the converse is not true: some graphs with zero star content have no articulation star for orders greater than 
5. An articulation star, otherwise known as a clique separator in the mathematical literature, is a completely 
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connected subset of vertices of a graph whose removal disconnects the graph. A biconnected graph is a graph 
with no articulation points, and so the set of graphs with no articulation star is a natural generalization. For 
the purpose of identifying graphs with an articulation star, we implement an algorithm [331 1341 which finds 
all articulation stars by finding a minimum ordering of the graph, and calculating the resulting fill-in graph. 
In practice, we used an algorithm which finds minimal rather than minimum orderings as this much faster; the 
trade off is that some graphs with an articulation star will not be identified, but for graphs generated in the 
calculation of Bg and Biq this occurred only about 2 percent of the time. Note that a graph is never falsely 
identified as having an articulation star when it does not. 

An alternative method would be to implement the algorithm of Whitesides in which a single articulation 
star is found in the same asymptotic time as all articulation stars are found by the minimal ordering method. 
However, the method of Whitesides perhaps should be explored as we are interested in small graphs, where the 
asymptotic behavior of the algorithm is perhaps not as important as the overhead involved in the calculation. 

Given the star content, and unlabeling factor, we may then proceed to calculate the virial coefficients up 
to order ten, in any integer dimension. The results of these computations have been given in Tabled 

4 Background 

We describe in this section the phase transition undergone by hard spheres in Subsection l4.1l and the rigorous 
information we have about this system in Subsection 14.21 Ideally, we would like to study the relation of the 
phase transitions to the virial coefficients. However the unambiguous interpretation of the ratio plots and 
their relation to the freezing transitions seen in computer studies is made difficult by the fact that we have 
very few exact results available for the hard sphere system beyond the general theorems that the pressure is 
continuous and non- increasing for positive densities in the physical region. For this reason in Subsection 14.31 
we discuss the various scenarios that have been proposed for the position of the leading singularity of hard 
spheres in order to provide a framework for later discussion in Section |5l 

4.1 Phase Transition 

Perhaps the property of hard spheres which is the most interesting and surprising characteristic is the existence 
of a phase transition, as discovered for D — 3 by Alder and Wainright A and Wood and Jacobson 5^ through 
molecular dynamics simulations in 1957, and for 13 = 2 in 1962 j^Tj. For dimensions three, four and five 
the phase transition is believed to be first order, while for hard discs the situation is more controversial despite 
intense research efforts over the past forty years |37j, j59j 4|1, 41 . Recently Jaster 42 showed that as the density 
of the system increases that hard discs first undergo a first or second order transition from the fluid phase to a 
hexatic phase with short range positional and quasi long range orientational order, and then undergo another 
second order transition to the solid phase which has quasi long range positional and orientational order. This 
is consistent with the Kosterlitz-Thouless-Halperin-Nelson- Young scenario, for which both transitions must 
be second order. 

The phase transition freezing {pf) and melting (ps) densities are given in Table 13 for D = 2 |39l I42|. 
D = 3 ^ESlEl! and — 4,5 |nHl- We fist also the density p^p and the scaled density B2Pcp = 2^~^rjcp of 
the densest lattices for dimensions D = 2, . . . , 8 from the lattice catalogue of Nebe and Sloane Note that 
the density "q — 1 corresponds to all space being filled. Finally, we include the lower bound of the radius of 
convergence obtained by Lebowitz and Penrose |46| . 
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Table 3: Values for i32, the density of the closest packed lattices, the densities at which freezing (p/) and 
melting (ps) occur, and the bound of Lebowitz and Penrose |46| . for hard spheres of diameter a in dimensions 
D = 2,--- ,8. 



D 


B2 


Pep 


B2Pcp — 2^ ^?7cp 


Pf/Pcp 


Ps/Pcp 


PLP/Pcp 


2 
3 


Tza^ 
2 

27rcr^ 
3 


2 

V3cr2 
%/2 


1.8137- •• 
2.9619- •• 


0.78 
0.66 


0.81 
0.75 


0.03990 
0.02444 


4 


7r (7 

4 


2 


4.9348 • • ■ 


0.50 


0.68 


0.01467 


5 


15 


2\/2 


7.4441 • • • 


0.41 


0.62 


0.00971 


6 
7 


7r (7 
12 

105 


8 

V3(t6 
_8_ 


11.9343- • • 
18.8990- ■ • 






0.00605 
0.00382 


8 


48 


16 

7^ 


32.4696- •• 






0.00224 



The existence of singularities foe I? > 3 at /?/ and ps is much more controversial than the singularities for 
D — 2. It has long been argued by Fisher 47^ that at these phase boundaries the pressure will be infinitely 
differentiable but will not be analytic and thus cannot be analytically continued into the region pf < p < ps 
from either side. These singularities have been proven to exist in the Ising model by Isakov |48) . but for 
hard spheres nothing rigorous has been proven. In this scenario the freezing density pf can be determined in 
principle from the leading singularity on the positive real axis of the low density equation of state. 

The alternative view of the first order (freezing) transition in hard spheres for _D > 3 is the assumption that 
there are no singularities at the phase boundaries and that analytic continuation from both sides is possible. 
The freezing transition is seen as a Maxwell construction making a convex envelope from a low and a separate 
high density free energy. In this scenario it is impossible to determine the freezing density from the low density 
of state alone. 

It is obviously of great theoretical importance to determine which of these two scenarios is correct for hard 
sphere for D > 3. This is a particularly difficult question if the radius of convergence of the virial series is 
determined by a singularity in the complex plane, which the ratio plots indicate is very likely the case for 
D > 4. In this case if the radius of convergence is less than the freezing density pf it is necessary to find a way 
to analytically continue the virial expansion beyond the radius of convergence to study a possible singularity at 
Pf. It is surely not possible to extract from our 10 term series both a leading singularity in the complex plane 
and to reliably continue the expansion beyond the radius of convergence to detect an infinitely differentiable 
singularity at pf. Consequently we will restrict our attention to locating the leading singularity in the complex 
plane. 



4.2 Rigorous Results 

There are few rigorous results available for the problem of hard spheres, and we will briefly summarize those 
few results here. 

Groeneveld jHHI proved that for the cluster expansion the radius of convergence must be greater than 
l/(2ei32), where e = 2.71828 • • • , and less than l/(2i?2). In addition, it is known that the cluster coefficients 
alternate in sign and hence the leading singularity is on the negative, real fugacity axis. Lebowitz and 
Penrose |46| adapted this bound to derive a corresponding expression for the virial series, given in Eq. 1231 

0.14476 



Plp > 



2B 



2 



0-14476 

VLP > -^jy- (23) 

Thus we know that the pressure is an analytic function of density for small densities; unfortunately in practice 
this bound does not seem to be close to the true radius of convergences as it is a long way from the density 
at which the phase transition is seen to occur, as can_^l^e seen from Table 13 



Later Fisher jHO] obtained bounds on the pressure in the vicinity of close packing for a D-dimensional 
system of hard particles. Extending the results of Hoover WD for hard parallel cubes, Fisher was able to prove 
that for T = ^ - 1 ^ that 

— < < — (24) 

where 0<<^i<l<^2<oo. For the general case of hard core particles, which includes hard spheres, he was 
able to establish the weaker relation 

for any 7?i < 1 < ?72- 



4.3 Scenarios for the Position of the Leading Singularity 

We present here various scenarios that exist for the location and nature of the leading singularity of the virial 
series for hard spheres in D dimensions. 

Over the past 40 years many approximate equations of state have been proposed for hard spheres, most 
commonly for the case D = 3. These approximates may be categorized by the location of their leading 
singularity. 

Many have high order poles at the space filling density rj — 1, including the solutions of the compressibility 
and virial Percus-Yevick integral equations j^j for £) = 3 by Thiele [22] and Wertheim |53l W>^ , the scaled 
particle theory of Reiss, Frisch, and Lebowitz [HI], a proposal by Guggenheim [HEI, and the widely used 
empirical formula of Carnahan and Starling |57j . 

The equations of state of Goldman and White and Hoste and Dael |^ have simple poles at or near 
the packing fraction of closest packed hard spheres. 

Other equations of state have a fractional power law divergence at or near the "random close packed" 
density r/rcp = 0.64 as defined by Ifill lf)3j . These approximates are obtained by constraining the 
divergence of the leading singularity to be of the form 

Pv/keT = A{tj - rjrcp)-' (26) 

As an example s is estimated as 1 in |H3 as 0.678 in and 0.76 in ^Hl- In |^ other values of rjrcp are 
chosen and the values of s lie in the range 0.6 < s < 0.9 depending on the approximation used. 

Some approximate equations of state include the freezing density, and we mention in particular that 
Torquato |68l I69| proposes an equation of state which agrees with the equation of Carnahan and Starling |57| 
for rj < rjf but which is of a different form for rj > rj-f. 

For dimensions greater than 3, we note the analytic solution by Leutheusser of the Percus-Yevick 
equation in D = 5, and the recent work of Robles, Lopez de Haro, and Santos |71| . in which the analytic 
solution was obtained for the Percus-Yevick equation in I? = 7. One of the most interesting results is that 
the radius of convergence is no longer determined by the singularity on the positive real axis at = 1, but 
instead by a singularity on the negative real axis for D > 5. 

Baram and Luban [HI in 1979 fitted the first seven virial coefficients using Levin approximants, and 
concluded that the leading singularity is at the close packed density. 

More recently, in 1994 Sanchez fitted the virial series with a Pade approximant, and then expressed 
the density in terms of the compression factor by inverting the series. He then fitted the inverted series with 
Fade approximants of increasing order to obtain estimates of the density at which the compression factor 
diverges. Sanchez then argued that as the order of the Pade approximants was increased that the position of 
the leading singularity converged to the density of close packing in two and three dimensions, and this was 
taken as evidence that the leading singularity of the virial series is in fact the density of closest packing. 

Gaunt and Joyce |74) provided a counter argument to the conclusion and methods of Baram and Luban, 
and indeed to Sanchez despite the fact that the work of Sanchez came well after. They performed a ratio 
analysis [JSI of the virial series, which they advocated as a robust method in the absence of any knowledge 
of the exact form of the function being described by the series. What was then the seven term virial series 



appeared to be behaving smoothly and ratio analysis suggested that the leading singularity lies on the positive 
real axis, perhaps at the space filling density. However, they pointed out that there is good reason to be 
cautious in extracting asymptotic behavior from only the low order coefficients: for several varieties of lattice 
gas with either nearest neighbor exclusion or first and second nearest neighbor exclusion, they show that the 
ratios initially behave smoothly indicating a singularity on the positive real axis, and then suddenly change 
behavior, and eventually change sign. This is because the asymptotic behavior is determined by a complex 
conjugate pair of singularities that mask the physical singularity on the real axis. 

The only models which demonstrate a phase transition for which exact results are available are hard core 
lattice gases in two dimensions. In particular the hard hexagon model [TSJ I77L 175] and a model of hard 
squares [79.1 have been exactly solved, and in both cases the radius of convergence is limited by a complex 
conjugate pair of singularities thus resulting in virial coefficients that oscillate in sign. 

Other relevant models include that of hard parallel cubes (HO! for which negative virial coefficients have 
been seen as low as sixth order, along with the Gaussian model |81[I82[E5| for which negative coefficients can 
be seen for dimensions D > 1. In particular for the Gaussian model oscillations are seen in the sign of the 
virial coefficients and this is evidence that the dominant singularities to are in the complex plane. 

Virial coefficients have been calculated for hard ellipsoids and other hard particles in two and three dimen- 
sions Uni- For these models we also see negative coefficients at relatively low order, with the effect being more 
dramatic for shapes which are far from spherically symmetric; in the case of ellipsoids, the effect is dramatic 
for large aspect ratio. 

The phase transition behavior of hard core systems ultimately depends on the geometry of the crystalline 
phase. Thus it is not immediately obvious how relevant hard hexagons and hard parallel cubes are to the 
hard sphere problem, as the crystalline phase fills all space and the phase transition has been observed to be 
second rather than first order. However, one could argue that this is only important in the vicinity of the 
phase transition itself, and that the leading singularities of the virial series for D > 2 may be a conjugate pair 
of singularities in the complex plane. 

5 Discussion 

The most important property of the virial coefficients Bk is not their actual numerical values for k less that 
some finite number but rather their asymptotic behavior as /c — > cxd because it is the asymptotic value which 
determines the radius of convergence. Of course no finite number of virial coefficients can give information 
on the k oo behavior unless there is some a priori reason to expect that the values of k are already in the 
asymptotic fc — > oo regime. 

In Subsection 15 . II we propose two such criteria for the asymptotic regime of hard spheres and will see that 
for fc < 10 there is no dimension in which both criteria are fully satisfied. Nevertheless it is still of interest to 
determine what results are obtained if well known methods are used in an attempt to determine the asymptotic 
behavior of the series from the first ten virial coefficients. In Subsections 15.21 and 15.31 we proceed with the 
methods of ratio analysis and of differential approximants respectively as robust and general methods for the 
purpose of series analysis. 

5.1 Criteria for asymptotic behavior of Bk 

Classical systems of hard particles for which the interaction potential is either infinite or zero, such as hard 
spheres and parallel hard cubes, are special in that some graphs have Ree-Hoover integrals which are identically 
zero for dimensions D less than some maximum value, as the restrictions imposed by the graph are so many 
that no configuration of points may satisfy them. Hence the number of Ree-Hoover integrals contributing at 
a particular order depends on dimension. We study this problem by performing the Monte-Carlo calculation 
outlined previously, but in addition we keep track of how many diagrams have been generated at least once 
from batch to batch whereas previously this information was thrown away. 

We tabulate our results for the number of contributing Ree-Hoover integrals in Table |21 For D = 1, only 
the complete star graph is non-zero, for D = 2 many graphs are zero at the orders calculated, but for D > 2 
it is difficult to make any estimate of how many zero diagrams there are. This is because there are a large 
number of graphs with numerically small Ree-Hoover integrals, and there is no way to estimate what the 



distribution of the values of these small diagrams is. Without this information it is impossible to estimate the 
number of zero diagrams until the Monte-Carlo procedure has run for a sufficiently long time such that almost 
all non-zero diagrams have been generated. It was found to be impossible to reach this regime for Bg and Biq 
in dimensions -D = 3, 4, and hence we suggest that new direct ways must be found to count zero diagrams if 
progress is to made on this problem for dimensions greater than two. 

The rate of growth of the number of non-zero Ree-Hoover integrals in two dimensions is far less than that 
of the number of biconnected graphs with non-zero star content. It is possible that the rate of growth has 
been reduced from 2'^('''~i)/^ to something like exponential growth, but there is not enough data to make a 
definitive statement on this issue. 

Criterion 1 

The number of nonzero Ree Hoover diagrams has approached its large k behavior. 

For = 10 this criterion is only completely fulfilled for D ~ 2 and is not fulfilled at all for D > 5. 
Our second criterion has been presented in our previous paper p2] 

Criterion 2 

The loose packed diagrams with the number of / bonds near their maximum value numerically dominate 

The validity of this criterion has been studied in detail in [22]. Here it was seen that for D = 3 and k > 12 
the criterion is well satisfied and that as D increases the criterion is satisfied for smaller values of k. However, 
for D = 2 the criterion was not satisfied even for k as large as 18. 

We thus conclude that there is no dimension in which both of these criteria are fully fulfilled although for 
D = 3 and Z) = 4 it is possible that they both could hold for some moderate values of k of the order of 12 to 
14. 

5.2 Ratio Analysis 

We begin our analysis of the virial coefficients given in Table by making a ratio analysis |74[ I75|. By 
comparison with the test function 

oo 

^akz'' with ak^kyz^, (27) 
fc=i 

which as 2; — > Zc has the singularity 

(l-z/z.)-l-^ (28) 

we see, up to possible logarithms which cannot possibly be seen in our ten term series, that the fc — > oo 
behavior of the ratios 

flfc+i/flfc - z-i(l + s/n + 0(n-2)) (29) 

indicates a simple pole at z = if s = and a divergence more (less) singular than a pole if s is positive 
(negative). Note in particular that if s < — 1 that the function is finite at z = and if s < — 1 — n the first n 
derivatives are also finite at z = Zj, even though there is a divergence in the (n -I- 1)"* derivative. 

We normalize the ratios to the density of the closest packed lattice by plotting BkPcp/ Bk-i and these 
ratios are plotted versus l//c in Figures l51ll2lfor D = 2, - ■ ■ ,8. The values of pcp obtained from the catalogue 
of lattice packings of Nebe and Sloane 45 are given in Table 01 

The qualitative description of our results is as follows. 

5.2.1 D = 2 

The plot of BkPcp/ Bk-i for = 2 is given in Fig. [S] The ratios are all positive and are smoothly decreasing. 
At fc = 10 the ratio is 0.9992 which is just less than the value of unity it would have if the leading singularity 
were at close packing. However, the ratios are still falling and hence suggest the radius of convergence is 
determined by a singularity at a density greater than close packing. For comparison the ratios for D = 3 are 
also plotted. 
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5.2.2 D > 5 



The ratios for D > 5 are plotted in Figures El and EH In these cases the first few virial coefficients are positive 
(giving positive ratios) and then to the order given the coefficients alternate in sign (giving negative ratios). 
This suggests that there is a singularity on the positive real axis which dominates the series initially, but there 
is a second singularity in the complex plane or negative real axis which is closer to the origin which determines 
the actual radius of convergence and dominates the ratios at high order. If the leading singularity is on the 
negative real axis the ratios will smoothly converge to some negative value, otherwise the ratios will oscillate. 



5.2.3 D^A 

The ratios for I? = 4 are plotted in Fig. ^2 There are no negative virial coefficients for fc < 10 but even 
though the accuracy for Bio is limited there seems to be a very definite oscillation which is developing just as 
oscillations develop for D > 5. Hence, from examining Fig. we expect that negative virial coefficients will 
occur for k not much greater than 12. 



5.2.4 D = 3 

The most important case is £> = 3 and this is also the most difficult to interpret. This case is plotted in Fig.El 
where it is compared with with D = 2 and is plotted on an expanded scale in Fig.^] In the plot of Fig.|51the 
ratios for D = 3 appear much the same as the ratios for D = 2 and extrapolate to a density significantly larger 
than Pep- However, an inspection of Fig. ^] reveals that unlike the situation for D ~ 2 the ratio plots for 
D = 3 are not always convex and can be seen from the magnitude of the slopes of the line segments given by 
5.8959 • • • , 4.43686(21), 2.80402(85), 2.212(32), 2.279(15), 1.745(75), and 1.34(34). This lack of monotonicity is 
evidence that the ratios for D — 3 are displaying small oscillations in amplitude. If these oscillations grow in 
magnitude then may be expected that eventually there will be oscillations in the signs of the virial coefficients 
and the radius if convergence will be determined by a singularity away from the positive real axis. 



5.3 Pade and Differential Approximants 
5.3.1 Notation 

The simplest way to make extrapolations of a finite number of terms in a power series expansion is the Fade 
approximant which assumes that the power series represents a rational function f(z) expressed as 

m - ^ (30) 

where Pm{z) and Qn{z) are polynomials of degree M and N respectively. We refer to this approximant as 
[M/N] . The coefficients of these polynomials are easily computed from the power series of M + — 1 terms 
by solving a system of linear equations This form cannot possibly represent a pressure which is continuous 
at its leading singularity and thus is not appropriate for investigating the question of equations singularities 
in equations of state for the fluid phase of any system. Nevertheless Fade approximates have been used by 
many authors |12[ 1141 115L [75| either to obtain information about the singularities of the virial series or as 
approximate equations of state for hard spheres. For reasons of comparison with the literature we list several 
Fade approximants for D = 2, 3 in the text below. 

A much more general method for extracting information from power series are the differential approximants 
presented in detail by Guttmann [7S]. In this method the power series is assumed to be represented by a 
function f{z) which is the solution of the linear differential equation 

Y,Q,{z;Li)z^ — f{z)^R{z-M). (31) 

1=0 

where Qi{z;Li) and R{z; M) are polynomials of degrees Li and M respectively 

Li M 

Q,{z■,L,)^Y.1k■.^z^ i?(0;M) ^^rfcz'^ (32) 

fe=0 IQ fc=0 



This form will incorporate the feature that P{p) /ksT = p + 0{p'^) for z — > by requiring that 



ro = 0, and ri ^ qo.i + (70,0 (33) 

The form H31|l is sufficiently general that it allows algebraic (or possibly logarithmic) singularities at the Lk 
zeroes Zj of Qk{z', Lk) The behavior near these zeroes f{z) in the case where the zeros are simple is 

/(z)^A(z,)|z-z,r- + ^(;-'f\ (34) 

with 

a,^K -I — — — (35) 

The special case K — 1 and R{z,M) = is called a Dlog Pade approximant and in this case f{z) is 
explicitly given as 

f(z) = zexp <^ - / dz \ (36) 

L Jo zQi{z;Li) J 

which near the zeroes Zi of Qi{z] Li) behaves as 

_ Oo('j:^o) 

f{z)^A{z,)\z-z,\ ^7«I(W. (37) 

This f{z) either diverges or vanishes at the singular points Zi and thus, like the Pade approximant the Dlog 
Pade is not able to accurately represent a singularity in the pressure at a point such as the freezing transition 
where the pressure is continuous. However, it should be adequate to study the question of whether for D — 
the leading singularity is in the complex plane off the positive real axis. 



5.3.2 Method 

We fitted the series with regular Pade, Dlog Pade and differential approximants, and attempted to get some 
measure of how well these approximants represented the series by seeing how effective they were at "predicting" 
coefficients when given a shortened series. This was done systematically by evaluating a measure which would 
quantify the relative mean square deviation of predicted coefficients from the actual coefficients of TablcQ] Wc 
reached the conclusion that there was no difference between Pade, Dlog Pade and differential approximants 
in extrapolating the virial series for hard spheres given the coefficients Bi, - ■ ■ , Bg. 

The Fortran program NEWGRQD given by Guttmann 75 was used to determine the differential approx- 
imants. We will now discuss the information that can be drawn from the approximants in each dimension in 
terms of the leading singularities, expected asymptotic behavior, and predicted coefficients. 

The scatter in the position and exponent of singularities, and in predicted coefficients comes from two 
sources: uncertainties in the virial coefficients, and differences between the approximants themselves. We 
take in to account the virial coefficient uncertainties by generating a reasonably large set (1000) of series 
by sampling from Gaussian distributions with the center given by the best values listed in Table ^ and the 
width given by the corresponding uncertainties. We then determine the Pade and differential approximants 
corresponding to each of these series, and combine the results to get mean values with uncertainties for the 
predicted coefficients of each approximant. In addition we obtain uncertainties for the coefficients of the 
rational function Pl{p)/Qm{p) for the Pade approximants. We found that the scatter in the predicted series 
coefficients due to uncertainties for Pade approximants matched the scatter among different approximants quite 
closely. The differential approximants were considerably more sensitive to the uncertainty of the coefficients, 
and generally were stable for around 2 to 3 less coefficients than the Pade approximants. Below we will 
neglect the discussion of scatter due to uncertainties in the coefficients and only consider the scatter between 
approximants. 

We list the singularities and corresponding exponents for the differential approximants for D — 2,3 in 
Tables m and [71 of Appendix IeI and the singularities and residues of the Pade approximants for D = 2,3 in 
Tables IHl and El of Appendix jFl Pade approximants are of little use in determining the singularities of a series 
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unless we have reason to expect the series to be weU represented by a rational function. As we have no reason 
to expect this to be the case for the virial series of hard spheres we do not make reference to these tables below, 
but include them for the interested reader to compare with the work of other authors |12[ll4l[T51l73| . We have 
not included tables for D = 4, • ■ • , 8 for reasons of length, and feel that the brief summary we make of the 
information in each of the tables is sufhcient. The dedicated reader may obtain the tables for singularities and 
exponents/residues in all dimensions, as well as tables of coefficients predicted by the approximants, either by 
downloading them |83| or upon request from the authors. 

Some care needs to be taken in interpreting the tables of singularities and exponents from the differential 
approximants, and we now discuss some of the relevant issues. Defective approximants, which have a singularity 
with an anomalously small exponent are marked with the symbol f , and are somewhat arbitrarily chosen as 
it is not possible to rigorously define what is an anomalously small exponent. We use the ad hoc procedure of 
Hunter and Baker [H3 j outlined by Guttmann [2^1 , that approximants be deemed defective if they have a pole 
with residue less than 0.003 that lies approximately inside the radius of convergence of the series. Defective 
approximants should be avoided as they are really lower order approximants in disguise, with an extra factor 
from a denominator pole and zero in the numerator that almost cancel. One other observation that should 
be made is that we expect in general if a series in the variable z has a singularity with exponent a, then 
the approximants may be expected to determine with more accuracy than a. Thus we expect more scatter 
between approximants in exponent values than in the position of singularities. In particular, even though we 
make an observation below about the position of the leading singularity in Tabled for D = 3, nothing can 
be said about the exponent. Further evidence for the unreliability of the exponents comes in the form of the 
large imaginary part of many of the exponents, which is very unlikely to be a true feature of the series. See 
Guttmann for a much deeper discussion on both of the issues raised above. 

We have combined the coefficients predicted by the various approximants [53| to produce Table 01 of 
predicted coefficients for dimensions D = 2, - ■ ■ ,8. We discuss these predicted coefficients as well as information 
obtained about singularities from the differential approximants for dimensions D = 2, - ■ ■ ,8 below. 



5.3.3 Analysis 

For D ~ 2, the leading singularities of the differential approximants lie on the positive real axis with the 
exceptions of [3,3;1] and [3, 2,2;0]. In most cases the singularity is located close to B2P — 1.98, with an 
exponent of —1.74(3). We list the singularities of the approximants in Table|51of Appendix IeI To make this 
estimate we have looked at all approximants which were not defective that use all 10 terms in the series, and 
calculated an average value. We discarded the values from the approximant [2, 2, 2; 0] which appeared to be 
an outlier. No other singularities follow a consistent pattern, and thus we see no signs of singularities that 
could eventually lead to sign changes in the virial series. As can be seen jSH], predicted coefficients are stable 
to Bis, and are listed in Table 01 where the scatter in the predicted coefficient is in the last digit. The Pade 
approximants [4/5] and [5/4] may be used as approximate equations of state at low density, and although they 
may be obtained directly from the data we include them for reference purposes along with the compilation of 
tables Pj . 

1 + 0.69939247(B2P) - 0.33033017(S2P)^ + 0.11294457(B2P)^ - 0.012320562(B2p)*' 
/[4/5l(P) - ^ _ o.30060752(B2P) - 0.81172709(B2p)2 + 0.62751627(B2P)^ - 0.17862580(B2P)'' + 0.021359218(B2P)5 



/[5/4] (P) ■ 



1 - 0.062894522(B2p) + 0.13851476{B2P)^ + 0.0067699403(B2P)^ + 0.0039942056(52 p)'* + 0.00047760798(B2P)^ 



1 - 1.0628945(B2P) + 0.41940485(B2P)^ - 0. 11367848(^2 p)^ + 0.021846467(B2P)'' 

(39) 

The variation in the polynomial coefficients in these approximants when one takes in to account the 
uncertainty in the virial coefficients is of the same order as the coefficients themselves. Hence one should 
not ascribe too much importance to the exact value of the Fade coefficients as they will change as future 
improvements are made in the accuracy of the virial coefficients. 

For D = 3, the position of the leading singularity varies more than for D = 2: in most cases the leading 
singularities for the high order differential approximants are in the complex plane with negative real part. 
However, the position and nature of this singularity varies substantially between approximants, and it is not 
possible to make a definitive statement as to whether this is a true singularity of the series, or an artifact of 



the approximant method. For approximants which only allow for 2 singular points, the leading singularity 
is on the real axis, and may be positive or negative. There does however appear in all approximants to be 
a stable singularity on the positive real axis, which is generally not the leading singularity. It is located at 
B2P = 3.75(3), with an exponent of —2.10(8), with the estimate made from all high order, non-defective 
differential approximants in Table [7| of Appendix IeI In accordance with the discussion of Subsection 15 . 21 this 
is the singularity that dominates the series initially. If the complex conjugate pair of singularities are indeed 
the leading singularities of the series then this will eventually lead to a change in sign of the virial series in 
3 dimensions. This can be confirmed or rejected by calculating more coefhcients, and possibly by increasing 
the accuracy of the coefficients to Biq which may reduce the scatter in the position of the singularities. As 
can be seen |H21, predicted coefficients are stable to Bie, and estimates of their values are given in Table 0] 
We explicitly give here the Pade approximants [4/5] and [5/4] which may be used as approximate equations 
of state at low density. 

1 + 0.50998996(B2p) + 0.20874890(B2P)^ + 0.036450422(B2P)^ + 0.0035176485(52 p)"* 
J[i/5](P) - 1 _ 0.49001003(62 p) + 0.073758936(B2p)2 - 0.018001749(S2p)3 + 0.0057761992(B2p)* - 0.00054759070(B2p)5 

(40) 

1 + 0.80408617(B2p) + 0.46662045(B2p)2 + 0.14197965(B2P)^ + 0.029738518(B2P)* + 0.0028192477(B2p)5 
/[5/4l(P)- 1 - 0.19591382(_B2P) + 0.037534279(B2P)^ - 0.060057986(^2 p)^ + 0.012302957(52 p)** 

As is the case for D = 2, if one takes in to account the uncertainty of the virial coefficients then the 
variance of the coefficients in the Pade approximants is of the same order as the coefficients themselves. 

In D = A there is no clear behavior for the leading singularity that can be gleaned from the approximants; 
there are cases where the singularity is on the positive real axis, on the negative real axis, or in the complex 
plane. Predicted coefficients are correspondingly variable, and this may be attributed to the low magnitude 
and correspondingly poor accuracy of Bg and especially -Bio- Almost all approximants predict that Bn will 
be positive (see Table ^J, but after that the predictions diverge: some allow for negative coefficients as soon 
as i?i2, and many allow for negative Bi^ or -B14, but this is in no way an improvement on our qualitatively 
based prediction in Subsection 15.21 that we will soon see negative coefficients for D = 4. Higher order virial 
coefficients and better accuracy for the known virial coefficients will help to resolve this question. 

For D — 5 the approximants give no clue as to the position or nature of the leading singularity. Predicted 
coefficients are stable only to B12, and are given in Table 01 Some approximants predict negative coeffi- 
cients for i?i3, and there are no predictions of asymptotic behavior that are consistent between the different 
approximants. 

We will discuss the cases of D = 6, 7, 8 together, as they are all very similar. The leading singularity is 
generally on the negative real axis, which suggests that the alternation in sign of the virial coefficients will 
continue for higher orders. There is no apparent pattern in the position of the next to leading singularity. 
Predicted coefficients are stable to i?i6 for D = 6, and to Bis for D = 7,8. This behavior is consistent with 
either the leading singularity being on the negative real axis as the differential approximant analysis suggests, 
or that there are complex conjugate singularities adjacent to the negative real axis but close enough so that 
there position cannot be accurately determined with ten virial coefficients. 



Table 4: Predicted coefficients for approximants with 10 exact coefficients in dimensions D — 2, ■ ■ ■ ,8, with 
uncertainty in the last digit. 



D 




-B12/-B2' 


^13 




Predicted coefficients 
Bi^/B},'^ Bi: 




Big 




-8177-82'^ 


Bis/bI' 


2 


1.089 X 10"^ 


5.90 X 10"^ 


3.18 : 


X 10"-' 


1.70 X 10"'^ 


9.10 


X 10"" 


4.84 : 


X 10"" 


2.56 X 10"" 


1.36 X 10"" 


3 


1.22 X lO"*' 


3.64 X 10"^ 


1.08 : 


X 10"^ 


3.2 X 10"** 


9.2 : 


X 10"^ 


2.6 > 


; 10"^ 






4 


1.2 X 10"" 






















5 


3 X lO^'^ 


-4 X 10"^ 




















6 


4.32 X 10"'' 


-3.68 X 10"" 


3.2 > 


: 10"" 


-2.9 X 10"" 


2.7 : 


X 10"" 


-2.5 


X 10"" 






7 


1.43 X 10"^ 


-1.40 X 10"^ 


1.42 : 


X 10"-'' 


-1.45 X 10"^ 


1.5 : 


X 10"^ 


-1.6 


X 10"^ 


1.8 X 10"^ 


-2.1 X 10"^ 


8 


2.56 X 10"^ 


-2.72 X 10"^ 


2.97 : 


K 10"^ 


-3.4 X 10"^ 


4.0 : 


X 10"^ 


-4.4 


X 10"^ 


5.5 X 10"^ 


-6 X 10"^ 
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6 Conclusion 



The key result is that negative virial coefhcients have been found for -D > 5, and a strong signal was found via 
the method of ratio analysis in Subsection IS . 2l that negative coefficients will soon occur for D ~ 4. On the basis 
of evidence from the methods of ratio analysis in Subsection 15 . 21 and differential approximants in Subsection 
15.31 there is a distinct possibility that at higher order some virial coefficients for = 3 will be negative. This 
would contradict the frequent assumption that all virial coefficients are positive for hard spheres in three 
dimensions. More work needs to be done to confirm or reject this hypothesis; in particular the calculation of 
Bii and possibly B12 may provide the necessary evidence. For D = 2 no deviation is seen from the behavior 
of a series that is dominated by a singularity on the positive axis at a density close to the space filling density. 
However, as discussed in Subsection 15.11 unless there is good reason to believe that the series has entered 
the asymptotic regime then it is impossible to draw strong conclusions about the position and nature of the 
leading singularity. 

Acknowledgments: This work was supported in part by the National Science Foundation under DMR- 
0302758. N. Clisby gratefully acknowledges support from the Australian Research Council. 

A Number of Configurations 

Numerical results for the virial coefficients are to be found in Tabled We list the number of batches of 10^ 
configurations in Table |S1 The total number of configurations can be obtained by multiplying the number of 
batches by the number of spanning trees used at a particular order. 

Table 5: Number of batches of 10^ configurations used in virial coefficient calculations, as a function of order 

and dimension 



D 


B4/BI 


B5/B2 


Be/Bl 


B7/BI 


Bs/Bl 


B9/BI 


Bio/Bl 


2 


1000 


9625 


9000 


9000 


15384 


19553 


6149 


3 


1000 


52573 


53463 


63751 


64675 


87609 


151349 


4 


1000 


8454 


9400 


10299 


21400 


31903 


38699 


5 


23199 


8436 


8618 


8597 


15607 


21042 


15398 


6 


1000 


8423 


8600 


8542 


5899 


6300 


1229 


7 


23010 


8213 


8600 


8500 


5898 


6300 


1300 


8 


1000 


8209 


8500 


8493 


5763 


6265 


1300 



B Star Content Results 

Here we plot histograms of the star content for labeled graphs contributing to Bq, Bj, Bg, Bg, and Biq. 
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Figure 2: Histogram plot of number of labeled graphs versus star content for Bq. 
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Figure 3: Histogram plot of number of labeled graphs versus star content for Br. 
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Figure 4: Histogram plot of number of labeled graphs versus star content for Bs 
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Figure 5: Histogram plot of number of labeled graphs versus star content for Bg. 
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Figure 6: Histogram plot of number of labeled graphs versus star content for Bio 

Spanning trees used in the Monte-Carlo Algorithm 
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Figure 7: Minimal sets of spanning trees for Bk, fc = 5, • • • ,10. 
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D Ratio Plots 
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Figure 12: Ratio plot for virial coefRcients in dimension D = 3 over a small domain to show non-monotonic 

behavior of the second derivative. 
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E Differential Approximants 



Table 6: Singularities for all differential approximants in £) = 2 in terms of B2P, with the corresponding 

exponents immediately below. Defective approximants are marked with j", and the singularities are arranged 
so that the singularity nearest B2P = 1.98 is in the left most column. In most cases, this singularity has the 
smallest modulus, but when this is not the case the singularity with smallest modulus is marked with *. 



Approximant 



Singularity / Exponent 



[3,2,2;0] 
[2,3,3;0] 
[2,2,2;1] 



[4, 3; 0] 


1.98 


1.58 ± 2.94 i 






-1.74 


0.528 =F 0.449 i 




[3, 4; 01 


1.98 


1.58 ± 2.92 i 


-30.6 




-1.74 


0.508 T 0.453 i 


0.728 


[4,4;0]1" 


1.99 


1.46 ± 2.80 i 


-1.11* 




-1.79 


0.406 =F 0.386 i 


3.18 X 


[5,4:01 + 


1.99 


1.46 ± 2.81 i 


-0.978* 




-1.79 


0.418 =F 0.390 i 


1.38 X 10"^ 


[4 S-Ol^ 


1.99 


1.46 ± 2.81 i 


-267. 




-1.79 


0.418 =F 0.390 i 


0.952 


\3 2- 11 


1.96 


—5.22 






-1.66 


12.3 




[2, 3: 11 


1.96 


-3.66 ± 7.70 i 






-1.67 


1.92 =F 7.72 i 




[3, 3: 11 

L") "-'5 -^J 


1.80 ± 0.198 i 


3.77 






-1.07 =F 0.0804 i 


-2.70 




[4 3- 11 


1.98 


-2.58 


-0.601* 




-1.74 


-6.19 


27.3 


[3, 4: 11 


1.98 


-0.589 ± 4.37 i 


-0.828* 




-1.76 


0.0841 =F 2.37 i 


8.79 


[2, 2: 21 


1.96 


-5.78 






—1.66 


14.0 




[3, 2; 2] 


1.95 


0.619* 






-1.59 


-0.935 




[2, 3; 2] 


2.06 


2.09 ± 0.978 i 






-2.25 


-0.830 ± 1.11 i 




[3, 3; 2] 


1.97 


0.945 ±1.80 i 






-1.70 


-0.344 ± 2.38 i 




[2, 2; 3] 


1.98 


0.345* 






-1.79 


-23.7 




[3, 2; 3] 


1.98 


0.361* 






-1.78 


-21.2 




[2, 3; 3] 


1.98 


42.6 


0.445* 




-1.78 


-3.23 


-16.6 


[2, 2; 4] 


1.98 


-0.0937* 






-1.75 


95.9 




[2,2,2;0] 


1.79 


1.21 






1.24 


-10.7 





2.16 ± 0.319 i 
-2.74 =F 0.955 i 
1.94 
-1.41 

1.98 
-1.72 



0.939 ± 1.61 i* 
-0.985 ± 1.12 i 
9.17 
-17.5 



-0.977* 
1.37 X 10"^ 
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Table 7: Singularities for all differential approximants in D = 3 in terms of B2P, with the corresponding 
exponents immediately below. Defective approximants are marked with f , and the singularities are listed from 
left to right in order of their modulus. The most stable singularity is on the positive real axis in the vicinity 
of B2P = 3.75, and this appears in the second column in all cases. 



Approximant 



Singularity / Exponent 



[4,3 


0] 


-1.03 ± 2.64 i 


3.71 






0.640 =F 0.0898 i 


-2.04 


[3,4 


0] 


-1.05 ± 2.73 i 


3.83 






0.752 =F 0.134 i 


-2.33 


[4,4 


0] 


-1.04 ± 2.65 i 


3.73 






0.652 =F 0.0885 i 


-2.09 


[5,4 


0] 


-1.04 ± 2.65 i 


3.73 






0.651 =F 0.0881 i 


-2.09 


[4,5 


0] 


-1.04 ± 2.65 i 


3.73 






0.651 =F 0.0881 i 


-2.09 


[3,2 


1] 


-3.49 


4.04 






6.71 


-2.95 


[2,3 


1] 


-1.68 ± 1.51 i 


3.79 






1.72 =F 0.928 i 


-2.25 


[3,3 


1] 


-1.93 ± 1.76 i 


3.81 






1.54 =F 1-36 i 


-2.30 


[4,3 


1] 


-1.11 ± 2.69 i 


3.73 






0.627 =F 0.190 i 


-2.07 


[3,4 


1] 


-0.621 ± 1.87 i 


3.79 






0.913 ± 1.18 i 


-2.24 


[2,2 


2] 


1.15 


3.64 






-3.22 


-1.99 


[3,2 


2] 


-0.544 


3.87 






13.4 


-2.47 


[2,3 


2] 


-1.24 ± 2.13 i 


3.80 






1.13 ± 0.230 i 


-2.27 


[3,3 


2] 


-0.464 ± 2.25 i 


3.78 






0.476 ± 1.21 i 


-2.22 


[2,2 


3] 


2.91 


3.75 






-1.17 


-1.96 


[3,2 


3] 


2.74 


3.77 






-1.44 


-2.01 


[2,3 


3] 


2.37 


3.77 






-2.05 


-2.09 


[2,2 


4] 


2.78 


3.77 






-1.38 


-2.00 


[2,2,2;0] 


4.35 


-5.19 






-3.86 


4.07 


[3,2,2;0] 


1.21 


6.04 






-1.35 


-9.45 


[2,3,3;0] 


-0.840 ± 2.45 i 


3.76 






0.469 ± 0.0986 i 


-2.17 


[2,2,2;1] 


-2.68 


3.70 






6.49 


-2.06 



-44.4 ± 
0.394 =F 



-6.75 
0.824 

-65.6 
14.1 

-232. 
183. 

58.5 i 
9.61 i 



-6.95 
-0.506 



-38.7 
-3.71 
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F Fade Approximants 



Table 8: Singularities for all Pade approximants in Z) = 2 in terms of B2P, with the corrcspoucliiig residues 
immediately below. Defective approximants are marked with f . 



Approximant Singularity / Residue 



[4/3] 


1.87 


2.40 


19.9 






-11.5 


13.7 


-63.2 




[3/4] 


1.91 


2.21 


2.78 ± 2.62 i 






-16.8 


17.4 


-0.568 ± 0.424 i 




[4/4] 


1.89 ± 0.187 i 


2.81 ± 1.43 i 








-1.33 T 6.26 i 


0.839 ± 1.90 i 






[5/4] 


1.90 


2.28 


0.513 ± 3.21 i 






-14.4 


15.8 


0.0116 T 0.0580 i 




[4/5]t 


-0.908 


1.94 


2.12 


2.60 ± 2.40 i 




-4.83 X 10"^ 


-26.8 


27.0 


-0.380 ± 0.551 i 



Table 9: Singularities for all Padc approximants in _D = 3 in terms of B2P, with the corresponding residues 
immediately below. Defective approximants are marked with f . 



Approximant 


Singularity / Residue 




[4/3]t 


0.516 


3.44 


3.92 




2.84 X IQ-"^ 


-146. 


178. 


[3/4] 


3.32 ± 0.602 i 


0.249 ± 7.66 i 






7.42 T 36.9 i 


-1.87 T 5.40 i 




[4/4] 


-2.65 


3.50 ± 0.461 i 


-16.3 




0.0100 


11.9 T 65.0 i 


-37.2 


[5/4] 


-1.13 ± 2.24 i 


3.57 ± 0.366 i 






0.00100 T 0.00669 i 


14.0 T 92.3 i 




[4/5] 


3.43 


-1.80 ± 3.99 i 


4.39 




-109. 


0.135 T 0.293 i 


213. 
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